from math import *
import pylab as pl
import random

def geo_scale(lat, lon):
	M = 6367450.
	phi = lon/360.*pi*2
	psi = lat/360.*pi*2
	lat_meters = 111133-560*cos(2*phi)+1.2*cos(4*phi)
	lon_meters = pi*M*cos(psi)/180.
	ratio = lat_meters/lon_meters
	return (lat_meters, lon_meters)

#call geo_scale first and then let
#let x = lat_meters*lat and y = lon_meters*lon

def total_curvature(x, y):

	n = len(x)

	total_k = 0
	length = (x[0]**2+y[0]**2)**0.5+(x[n-1]**2+y[n-1]**2)**0.5

	for i in range(1, n-1):

		t1x = x[i+1]-x[i]
		t1y = y[i+1]-y[i]

		t0x = x[i]-x[i-1]
		t0y = y[i]-y[i-1]

		s1 = (t1x**2+t1y**2)**0.5
		s0 = (t0x**2+t0y**2)**0.5

		T1x = t1x/s1
		T1y = t1y/s1

		T0x = t0x/s0
		T0y = t0y/s0

		s = 0.5 * ((x[i+1]-x[i-1])**2+(y[i+1]-y[i-1])**2)**0.5

		kx = (T1x-T0x)/s
		ky = (T1y-T0y)/s

		k = (kx**2+ky**2)**0.5*s

		total_k += k
	return total_k

def length(x,y):
	n = len(x)

	length = 0

	for i in range(0, n):
		length += (x[i]**2+y[i]**2)**0.5



	
'''
#test code

x = []
y = []
z = []

N = 100
R = 10.
h = 10
t = []
s = 0

for i in range(0, N):
	u = random.randint(1, 10)
	t.append(s+u)
	s += u 

for i in range(0, N):
	t[i] = t[i]*1.0/s


for i in range(0, N):
	x.append(R*cos(1.*t[i]*2*pi))
	y.append(R*sin(1.*t[i]*2*pi))
	z.append(t[i]*h*2*pi)

print total_curvature2(x,y,z)-R/(R**2+h**2) * (R**2+h**2)**0.5 *2*pi
'''
	

